Tuning MOF/polymer interfacial pore geometry in mixed matrix membrane for upgrading CO2 separation performance

The current paradigm considers the control of the MOF/polymer interface mostly for achieving a good compatibility between the two components to ensure the fabrication of continuous mixed-matrix metal-organic framework (MMMOF) membranes. Here, we unravel that the interfacial pore shape nanostructure plays a key role for an optimum molecular transport. The prototypical ultrasmall pore AlFFIVE-1-Ni MOF was assembled with the polymer PIM-1 to design a composite with gradually expanding pore from the MOF entrance to the MOF/polymer interfacial region. Concentration gradient–driven molecular dynamics simulations demonstrated that this pore nanostructuring enables an optimum guided path for the gas molecules at the MOF/polymer interface that decisively leads to an acceleration of the molecular transport all along the MMMOF membrane. This numerical prediction resulted in the successful fabrication of a [001]-oriented nanosheets AlFFIVE-1-Ni/PIM-1 MMMOF membrane exhibiting an excellent CO2 permeability, better than many MMMs, and ideally associated with a sufficiently high CO2/CH4 selectivity that makes this membrane very promising for natural gas/biogas purification.


Supplementary Text
The initial step involved the optimization of the crystal structure of AlFFIVE-1-Ni by using Density Functional Theory (DFT) calculations, allowing for complete relaxation of both the atomic positions and cell parameters.All these calculations were performed using the Vienna ab initio Simulation Package (VASP, version 5.4.4) (40) with a plane wave energy cut-off of 650 eV and convergence criteria of 10 -5 eV and 10 -2 eV/Å for the energy and sum of the atomic forces, respectively.The electronic exchange-correlation interaction was treated using the Perdew-Burke-Ernzerhof (PBE) functional within the generalized gradient approximation (GGA).(39) The first Brillouin zone was set with a 3×3×2 Monkhorst-Pack k-point mesh.The van der Waals correction DFT-D3 method (53) was also used to accurately account for long-range dispersion interactions.In order to better describe the d states of the Ni of AlFFIVE-1-Ni, an Hubbard U correction of 6.4 eV was used.(54) Next, the DFT-optimized crystal structure of AlFFIVE-1-Ni model was further employed to construct the AlFFIVE-1-Ni (001) surface model.The plane (001) was chosen to orient the cleavage of the surface to mimic the experimental scenario corresponding to the fabrication of AlFFIVE-1-Ni nanosheets with exposed (001) facets experimentally reported earlier.(21) The slab model has a zero dipole in the direction perpendicular to the surface slab terminated by F-atoms coordinated to the Ni-atoms, as shown in Fig. S1.This surface model was then geometry-optimized at the DFT level using the same level of theory and parameters as for the optimization of the bulk model.The density derived electrostatic and chemical (DDEC06) method, as implemented in the CHARGEMOL module,(55) was employed to calculate the partial atomic charges for the surface model.
A final surface slab model was then constructed for the force field simulations by extending 5 times the surface model in the x and y directions giving final x and y lengths of 48.925 Å.In these computations, the surface slab model was regarded as a fully flexible framework employing intramolecular potential parameters taken from the Universal Force Field (UFF) (43).Intermolecular interactions were modeled using 12-6 Lennard-Jones (LJ) van der Waals and Coulombic contributions.The cutoff is set to 12 Å for all following simulations.The 12-6 LJ parameters were taken from the generic forcefields UFF (47) and DREIDING (44) for the atoms of the inorganic and organic nodes respectively.The 6FDA-DAM polymer model utilized in this study was adopted from our previous work.(56) Initially, it was constructed by starting with a united-atom monomer, as depicted in Fig. S2, and undergoing a controlled polymerization process (61).A comprehensive description of this methodology can be found in our earlier publications.(28,57), The terminations of the resulting polymer were achieved as follows: on one end, the replacement of a nitrogen atom by an oxygen atom and, on the other end, the addition of -NH2 group to a carbon atom (Fig. S3).This polymer model was treated as fully flexible in our calculations: bonded interactions were modeled using the analytical expression and parameters described by the generic forcefield GAFF while the nonbonded interactions were treated by a LJ potential contribution and a coulombic term, with parameters taken from our previous work and reminded in Table S1.S1 12-6 LJ parameters and charges for 6FDA-DAM model.
To construct the PIM-1 polymer model, the Polymatic code (57) was also utilized in conjunction with the LAMMPS code (41,42).The PIM-1 configuration was derived from this initial model by performing annealing cycles via MD simulations in two NVT ensembles: NVT high and NVT low .
T high and T low were varied between 100 K and 900 K to achieve diverse equilibrium configurations.
Atom types, LJ parameters and charges for the PIM-1 model.Supplementary Fig. S3) and Nitrogen (nitrile group, cf.The gradient was maintained with adaptive forces located in IFR (inlet force region) and OFR (outlet force region) that enable to control the concentration inside the control volumes (ICR and OCR) to keep a target value.This is achieved by imposing two different target values in the inlet and outlet control region, where the forces are defined by the following equations:

Table S3
The parameters used in the permeation simulations.ZF is the centre of external biasing force located region, w is the width of the external biasing force region and ki is the force constant.S5 Exponential fit data for rotational autocorrelation functions of CO2 at the heterojunction of both composites reported in Figure 3d.The rotational autocorrelation functions are calculated via C2(t) = <P2(e(t).e(0))>where P2 is Legendre polynomial of second-rank and e(t) is the unit vector corresponding to orientation of the molecules.The reorientation time scale of CO2 in both cases can be obtained by fitting C2(t) by an exponential decay function (  -.% ).The corresponding parameters as well as the fitting error R are reported.Relevant time scale for librational motion corresponds to 1/A.S7 Evaluation of the single-file mobility factor F at difference z-distance of the MOF for ALFFIVE-1-Ni embedded in the AlFFIVE-1-Ni/PIM-1 composite model using a lag-time approximation adapted to single-file diffusion  Single gas feed pressure: 2 bar; permeation at 35 °C. 1 Barrer = 10 -10 cm 3 (STP) cm/cm 2 s cmHg.

Table S8
A comparison of the gas separation performances and permeability (P) enhancement of

Fig. S24
A plot of the CO2/CH4 selectivity vs CO2 permeability for the different MOFs/PIM-1 polymerbased MMMOFs reported in the literature (references are provided in Table S10).
Fig. S1 provides an illustration of the atom types present in the AlFFIVE-1-Ni-(001) surface model.(cf.ref.21 for the respective charges and non-bonded parameters).

Fig. S1
Fig. S1 Atom types associated with the AlFFIVE-1-Ni (001) surface model viewed from (a) side and (b) top, respectively.Only a single end of the surface slab structure is shown on the left for clarity.F, N, C, Ni, Al and H are respectively shown in green, blue, gray, orange, purple and white, respectively.This figure highlights that the MOF surface is terminated by F-atoms.

Fig. S2 Optimized
Fig. S2Optimized 6FDA-DAM monomer structure, with head and tail atoms represented by atoms N17 and C30.

AlFFIVE- 1 -
Ni/6FDA-DAM and AlFFIVE-1-Ni/PIM-1 composite models were obtained by using the computational strategy we previously developed and validated on a series of MOF/Polymer systems.First, the atomic coordinates of the polymer model underwent unwrapping in the z-direction, and subsequently, the simulation box was adapted to conform to the lattice parameters of the AlFFIVE-1-Ni (001) structure model, which measured 48.925×48.925×150.000 in dimensions.Subsequently, the polymer model was brought into contact with the AlFFIVE-1-Ni (001) slab model, and subjected to a 21-step molecular dynamics (MD) simulation regimen, which encompassed NVT and NPT cycles.These cycles were controlled using the Berendsen thermostat and barostat, with relaxation times of 0.1 ps and 0.5 ps respectively, within a customized DL_POLY Classic code.(59,60)This method facilitated the polymer's equilibration, leading to the formation of a compact AlFFIVE-1-Ni/polymer composite.The interactions between the MOF surface and the polymer were treated by incorporating van der Waals and Coulombic terms, with the cross 12-6 LJ parameters determined using the Lorentz Berthelot mixing rule.(61)The two components of the composite were treated as fully flexible with the force field parameters mentioned in the previous sections.The van der Waals interactions were handled with a 12 Å cutoff while the electrostatic interactions were calculated using the Ewald summation method with a tolerance of 10 -6 .(62)These calculations were performed during a total simulation time of 10 nanoseconds (ns) with a time step of 1 femtosecond (fs).

Fig. S6
Fig. S6Calculated pore size distributions (PSD) for the two different composite models.The plot is statistically determined for each composite system five independent MOF/Polymer interface configurations were considered.The PSD calculation is performed for the composite excluding the contribution of the AlFFIVE-1-Ni.

Fig. S7
Fig. S7 Illustration of the free porosity over the xy-plane for AlFFIVE-1-Ni/6FDA-DAM composite model represented in Fig.1 for different z-values when one moves away from the MOF surface, where gray regions represent void spaces, and white regions denote the presence of polymers.

Fig. S8
Fig. S8 Illustration of the free porosity over the xy-plane for 5 different configurations of AlFFIVE-1-Ni/6FDA-DAM for different z-values when one moves away from the MOF surface, where gray regions represent void spaces, and white regions denote the presence of polymers.

Fig. S9
Fig. S9 Illustration of the free porosity over the xy-plane for AlFFIVE-1-Ni/PIM-1 composite model represented in Fig.1 for different z-values when one moves away from the MOF surface, where blue regions represent void spaces, and white regions denote the presence of polymers.

Fig. S10
Fig. S10 Illustration of the free porosity over the xy-plane for 5 different configurations of AlFFIVE-1-Ni/PIM-1 for different z-values when one moves away from the MOF surface, where gray blue regions represent void spaces, and white regions denote the presence of polymers.

Fig. S11 (
Fig. S11 (a-b) Radial distribution function calculated between CO2 and the Oxygen (carbonyl group, cf.Supplementary Fig. S2) and Nitrogen atoms (cf.Fig. S2) of 6 FDA-DAM at the interface of the AlFFIVE-1-Ni/6FDA-DAM composite.Simulations were carried out at 300K with 1 bar pressure in the corresponding composite.

Fig. S12 (
Fig. S12 (a-b) Radial distribution function calculated between CO2 and the Oxygen (ether group, cf. Fig. S3) atoms of PIM-1 at the interface of the AlFFIVE-1-Ni/PIM-1 composite.Simulations were carried out at 300K with 1 bar pressure in the corresponding composite.

Fig. S13
Fig. S13 Atomic density plot of MOF and polymer as a function of z-coordinate for the representative configuration of AlFFIVE-1-Ni/6FDA-DAM composite.We remind that the interfacial region represented in dashed lines is defined by the limit between the z value for which the polymer atomic density starts to oscillate around an equilibrium state, and the z value of the center position of first metal(Ni) -pyrazine square-grid layer of AlFFIVE-1-Ni in the composite vanishes.

Fig. S14 Atomic
Fig. S14 Atomic density plot of MOF and polymer as a function of z-coordinate for the representative configuration of AlFFIVE-1-Ni/PIM-1 composite.We remind that the interfacial region represented in dashed lines is defined by the limit between the z value for which the polymer atomic density starts to oscillate around an equilibrium value and the z value of the center position of first metal(Ni) -pyrazine square-grid layer of AlFFIVE-1-Ni in the composite vanishes.

Fig. S15 Schematic
Fig. S15 Schematic representation of the simulation box used in the CGD-MD simulations.IFR=Inlet Force Region, ICR=Inlet Control Region, ITR=Inlet Transition Region; OTR=Outlet Transition Region, OCR=Outlet Control Region and OFR=Outlet Force Region).

Fig. S16
Fig. S16 Averaged pressure control at the (a) inlet and outlet (b) control regions for AlFFIVE-1-Ni/6FA-DAM composite model along the 500 ns CGD-MD simulations.

Fig. S22 MSD
Fig. S22MSD plot for CO2 obtained from EMD simulations performed at 1200 K and 10 bar CO2 pressure for 200 ns (100 ns equilibration followed by 100 ns of production) on the pristine MOF AlFFIVE-1-Ni averaged over 5 different MOF configurations.
Evaluation of the single-file mobility factor F at difference z-distance of the MOF for AlFFIVE-